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Multiparticle symmetrization effects are contributions to 
the spectra of Bose-symmetrized states which are not the 
product of pairwise correlations. Usually they are neglected 
in particle interferometric calculations which aim at deter- 
lining the geometry of the boson emitting source from the 
neasured momentum distributions. Based on a method in- 
roduced by Zajc and Pratt, we give a calculation of all mul- 
tiparticle symmetrization effects to the one- and two-particle 
momentum spectra for a Gaussian phase space distribution 
jjj emission points. Our starting point is an ensemble of N- 
fc-phrticle Bose-symmetrized wavefunctions with specified phase 
^space localization. In scenarios typical for relativistic heavy 
■^jen collisions, multiparticle effects steepen the slope of the 
,-uane-particle spectrum for realistic particle phase space densi- 
ties by up to 20 MeV, and they broaden the relative momen- 
C*F)|im dependence of the two-particle correlations. We discuss 
£$hese modifications and their consequences in quantitative de- 
tail. Also, we explain how multiparticle effects modify the 
S^rormalization of the two-particle correlator. The resulting 
^normalization conserves event probabilities, which is not the 
^jise for the commonly used pair approximation. Finally, we 
(~j»ropose a new method of calculating Bose-Einstein weights 
(«-*rom the output of event generators, taking multiparticle cor- 
relations into account. 
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^ACS numbers: 25.75.+r, 07.60.1y, 52.60.+h 
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I. INTRODUCTION 

> ■ 

Most hadrons are emitted in the final stage of a rela- 
Jtivistic heavy ion collision. They do not probe directly 
tne hot and dense intermediate stages where quarks and 
gluons are expected to be the relevant physical degrees of 
freedom for equilibration processes. The geometrical size 
and dynamical state of the hadronic phase space emission 
region, however, depends sensitively on the entire evolu- 
tion of the collision. This motivates current attempts 
to reconstruct its spatial and dynamical state from the 
experimental hadron spectra and to use it as a starting 
point for a dynamical back extrapolation into the hot 
and dense intermediate stages Two-particle corre- 
lations of identical particles which arc sensitive to the 
space-time characteristics of the collision [Q , play a cru- 
cial role in this approach. The reconstruction program 
based on their analysis has very good prospects: due to 
the increasing event multiplicities and larger statistics of 
the CERN SPS lead beam program (and the yet better 
quality data expected from Relativistic Heavy Ion Col- 



lider RHIC at BNL), particle interferometric measure- 
ments start showing statistical errors on the few percent 
level. Also, systematic errors are increasingly better un- 
derstood. Theoretical calculations should aim for a simi- 
lar accuracy and control necessary approximations quan- 
titatively. 

One uncontrolled approximation used so far in almost 
all particle interferometric calculations is to neglect for 
the particle momentum spectra of Bose-Einstein sym- 
metrized TV-particle states all multiparticle correlations 
which cannot be written in terms of simpler pairwise 
ones. This reduces the number of AH terms contributing 
to the two-particle correlator C(K,q) to a manageable 
sum over all particle pairs. Beyond this approximation, 
two approaches have been used in the literature. First, 
Zajc has employed Monte Carlo techniques [|[ to gener- 
ate events with realistic multiparticle correlations. This 
amounts to a shifting prescription {p^} — > {p^} which 
modifies the momentum distribution of simulated events 
according to unit weights (which themselves depend on 
the space-time structure of the source). Second, Zajc has 
found H the first steps towards a calculational scheme 
brought into final form by Pratt |4j : for model distribu- 
tions, the A^-particle spectra are given by a simple algo- 
rithm involving only two types of terms: C m and G m . In 
practice, this reduces the sums over all AM permutations, 
typical for the calculation of momentum spectra, to sums 
over all partitions of N. 

For both these approaches, there are first numerical 
calculations H-f| and related analytical attempts H [nj 
to control multiparticle effects to the one- and two- 
particle spectra, but a detailed study of their momentum 
dependence is missing, even for simple models. This work 
aims at filling this gap, making quantitative statements 
about the extent to which the slope of the one-particle 
spectrum and the relative momentum dependence of the 
two-particle correlations are modified due to multiparti- 
cle symmetrization effects. Our investigation takes the 
set £ of phase space emission points (pi,rj,ij) as initial 
condition. For notational simplicity, we restrict the dis- 
cussion to one particle species, like-charge pions say. To 
the set E, we associate a symmetrized A^-particle wave 



function 11 



' s£Sn \i=l / 

where the sum runs over all permutations s € <Sjv of 



1 



the N indices, X is a shorthand for the N 3-dimensional 
coordinates X,, and the functions fi denote single par- 
ticle wavepackets centered around (pi,^) at initial time 
U and propagated according to the free time-evolution. 
Final state interactions, which imply a structure of the 



TV-particle state different from (1.1), will not be consid- 
ered in the present work. The wave function \I/ n defines 
the boson emitting source. What we are interested in is 
the calculation of the one- and two-particle momentum 



spectra resulting from (1.1), the information they con- 
tain about the initial distribution of the 'emission points' 
z i = (Pi; r ii U), the extent to which these results modify 
the predictions of the pair approximation, and finally the 
algorithm which implements the numerical calculation of 
multiparticle spectra from the initial distribution E. 

Our work is organized as follows: Section |l| shortly 
sets up and illustrates the general formalism via which 
particle momentum spectra are calculated from an N- 
particle state. In section [II, we discuss the properties 
of Gaussian wavepackets w hich we choose for the single- 
particle states fi(Ki,t) in (1.1). Section IV contains the 
main results. It gives a complete qualitative and quanti- 
tative analysis of multiparticle contributions to the one- 
and two-particle spectra for a boson emitting source of 
Gaussian phase space distribution. In section we dis- 
cuss the implications of this model study for an algorithm 
which calculates the Hanbury Brown Twiss (HBT) ra- 
dius parameters and momentum spectra for an arbitrary 
distribution of emission points. Results and related con- 
ceptual questions are summarized and discussed in the 
Conclusions. 



II. THE FORMALISM 

We want to determine for the iV-particle symmetrized 
wave function 4"^ the detection probability for measur- 
ing N identical bosons at time t at the positions Xi and 
momenta P^. We calculate first the A^-particle Wigncr 
phase space density for Wjv IM 

3iV, 



W N (X,P,t) = (2n) iN * N (X,t) 



f N \ 

n<5 (3) (p*-pi) net*) 

VJ=1 / 



=m e n^^-pi.*). ( 2 - ia ) 

s,s'es„ i=i 

W l3 (X,P,t) = (2Tr) 3 f l (X,t)5^(P-P)f*(X,t) 

= |d 3 yc- tPy /,:(X+|,t)/;(X-f,i). (2.1b) 

Here, P denotes the momentum operator acting on 'J at. 
The one-particle pseudo-Wigner functions Wij(X,P,t) 
provide the basic building blocks for the calculation of 



the A^-particle momentum spectrum which is obtained 
by integrating (2.1a) over all spatial coordinates, 



Vn(P)=M* / d 3 XW N (X,P,t) 



N 



= W E (2.2a) 

s ,s'es„ i=i 

T l3 {P) = J d 3 XW l3 (X,P) = D t (P)D*(P), (2.2b) 

where A* is a normalization constant, A/* = 
1/(^n\^n), ensuring that the probability of detecting 
N particles is one. For a free time evolution of ^n, the 
integration over the spatial components of (2.2a) leads 
to a time- independent expression 'Pjv(P), since interac- 
tions between the particles are necessary to change the 
momentum distribution in time. In contrast, integrat- 
ing WV(X, P,t) over all momenta leads to the detec- 
tion probability of the N bosons at positions Xj , which 
is a time-dependent quantity since free evolving bosons 
change their positions in time. The calculation of the 
AT-particle momentum spectrum according to ( [2. 2a ) in- 
volves a sum over (AH) 2 terms. Due to the factorization 
of ^ (P), this reduces to a sum over AH terms, 



Vn(P) = 



A/* 
A! 



AT 



sES„ 1=1 



(2.3) 



In what follows, we are especially interested in the 
one- and two-particle momentum spectra 'P)y(Pi), 
■p^r (Pi,P2) associated with the A^-particle state ^jv- 
These are obtained by integrating Vn(P) over all but 
one, respectively two momenta, 



A' 



JV 



a; 



(2.4a) 



1=2 



N 



P&(Pi,P a ) = ^^^ lSl (Pi)^ S2 (P 2 )n^ S! > (2.4b) 



fa = J d 3 PT tJ {P). 



1=3 



(2.4c) 



All particle momentum spectra are given in terms of the 
building blocks Di(P) (which determine Tij) and fij. 
Once the analytical form of the single particle wavefunc- 
tions fi is specified, these are readily calculated. In what 
follows, capital letters denote measurable position and 
momentum coordinates, small letters p^, r i; tj denote 
the centers of wavepackets which are not directly mea- 
surable. The only exception to this is the measurable 
relative momentum q = Pi — P 2 of the two-particle cor- 
relator C(K = | (Pi +Pa)j q) which we denote by a small 
letter. 
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A. The Zajc-Pratt algorithm 

Dynamical correlations between particles in the source 
are reflected in correlations in the set of emission points 
(pi,rj,tj). If there are no correlations, then the initial 
distribution of the centers of single particle wavepack- 
ets is given by a one-particle probability distribution 
p(p, r, t). The n-particle spectra for a set of events with 
multiplicity N are obtained by averaging over this distri- 
bution 



Piv(Pi,..,P«) 

Vpi = d 3 Ti d 3 pi dti p(p u r t ,ti) . 



(2.5b) 



A particular model distribution with correlations which 
assumes that the emission probability of a boson is in- 
creased if there is another emission in its vicinity ]10| ], is 
e.g. obtained by replacing in (2.5a) 



A' 



JY 



n p(vt> r *' *<) — y n p ( pi > r *> **) 



w{N) 



(2.6) 



Here lu(N) is an averaged normalization defi ned be- 
low. The technical advantage of adopting (2.6) is that 
the (pj, rj, ij)-dependent normalization factor A/* in the 
spectrum of ( 2.5a ) is canceled. This allows to write 
without approximation all spectra in terms of the build- 
ing blocks 



g w (Pi,p 2 ) - J A'^PO/ii*. h 

x/ im - 1 . m A m (P 2 ), 
Jd 3 PG m (P,P). 



x .. 



(2.7a) 
(2.7b) 



The resulting Zajc-Pratt (ZP) algorithm for the calcula- 
tion of one- and two-particle spectra reads P,p|jl0| 



>w= E 



AH 



c^ct? ...c\ 



(n,l n )N 
N 



n n n'-(«n!) 



'In 
n J 



(2. 



v>) = E 



(N-l)\ w(N-m) 



^(N-m)\ w(N) 



G m (P,P), (2.9) 



A 



V 2 ( P p ) y r (N-2)\w(N-J) 

p w (P!,p 2 ) - w(n) 

.7-1 

xE [G,(Pi,Pi)Gfj_i(P 3 ,P 3 ) 



+G i (Pi,P 2 )G J _ i (P 2 ,Pi) 



(2.10) 



These spectra are normalized to unity. In Appendix |A|, 
we give their derivation in some combinatorical detail. 



Due to the ZP-algorithm, the calculation of the N- 
particle spectra is reduced from sums over all permu- 
tations to sums over all partitions (n 7 l n ) N of a set of 
N points into l n subsets of n points, N = J £2 n l n n. The 
number of partitions of N grows asymptotically like e*~^ ' . 
For explicit calculations with event multiplicities in the 
hundreds, it is hence important to get control over the 
m-dependence of the Pratt terms C m and G m . This is 



our strategy in section IV 



B. A simple example: the Zajc model 

To illustrate the above formalism, we consider a nor- 
malized iV-particle density matrix pi for multiparticle 
states |xi, ...,xjv), created by repeated operation of the 
single-particle creation operator $ (x) , 



^t(x) 



d 3 k 



■e lkx ff (k)4, 



(2tt) 3 /2' 

pi N) =N J dx 1 ...dx N p{x 1 )...p{x N ) 
x |xi, ...,xjv)(xi, ...,Xjv| ■ 



(2.11a) 



(2.11b) 



This density matrix specifies in particular the onc- 
and two-particle spectra Trfp^apflp] and Tt[p^ ap i a t P2 
apjQpJ. For the Gaussian model distribution p] 



l.g(k)| 2 
(2tt) 3 

P(x) 



(2^ ) -3/2 exp[ _ k 2 /2p 2 ] 

(^R 2 )- 3 / 2 exp[-x 2 /i? 2 ] , 



(2.12a) 
(2.12b) 



the effects of multiparticle correlations on the HBT- 
radius parameters have been considered already by Zajc. 
His discussion however is restricted to an explicit cal- 
culation of three-particle symmetrization effects and to 
qualitative estimates of higher order contributions. Here, 
we demonstrate that the ZP-algorithm allows for a com- 
plete quantitative analysis of the Zajc model. The first 
step is to identify the building blocks of the particle spec- 
tra (O), 



A(P) 



5(P) 



(2tt) 3 /< 



exp [iP ■ Xj] , 



fij = exp 



-f(x l -x J ) 2 



(2.13a) 
(2.13b) 



The calculation of the terms C m reduces then to Tri- 
dimensional Gaussian integrations. One finds G\ = 1 
and 



c m = (i+ P 2 i? 2 r 3 ( m - 1 )/ 2 



/ 1 -h 
\ det B rn 



3/2 



det B„ 



h = 



2m-i 
1 

1 + 1/c 



(Tm (}) 1) 



(2.14a) 



(2.14b) 
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where T m denotes Chebysheff polynomials of the first 
kind pq|. The momentum-dependent terms read 



G m (K,q) 



(m) 



(m) 

9k 



(1 + g ( K ] f/ 2 r K 2 (m) 21 
(27rp2)3/2 exp ^ 9 k / z Po\ 

x exp[- (i? 2 4" l) /4 + 1/8^) q 2 ] (2.15a) 

e+ , = + S im , (2.15b) 

1, eL 4) = <Jii - <5 lm , (2.15c) 
(2.15d) 



Is* /l-i 



(to) 



(A m L = (1 + c)<5y - § (6i, j+1 + 5 i+hj ) 



The main message of these involved but explicit expres- 

(r, 

sions is contained in the m-dependence of the terms gg 

and g^- These specify the K- and q-dependence of 
the building blocks G m and hence of all spectra. Here, 
they are functions of the phase space volume V — Pq i? 3 
only, and their behaviour can be understood by simple 
arguments: 
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FIG. 1. The terms Sq"', and the remainder term 



(1 — K)l det B m characterize the q-dependence, K- dependence 
and weight of Pratt terms and hence the momentum depen- 
dence of particle spectra. Their m-dependence contains infor- 
mation about how higher order multiparticle correla tions affect 
the spectra. Results are shown for the Zajc model (2.1i). 



The factors 



(m) 
A' 



generically increase with increasing 



to, see Figure yj. Especially, in a limiting case, one finds 

(2.16) 



hm g y K 

V—>oo 



m . 



The reason is that Bose-Einstein symmetrization effects 
enhance the low momentum region of the one-particle 
spectrum, leading to steeper slopes. This one-particle 
spectrum V N (P) is a linear superposition of Gaussian 
terms G m (K, q = 0), which due to g^ show increasingly 
steeper slopes. 

The q-dependence of G m governs the q-dependence of 
the two-particle spectrum. The terms g^ depicted in 
Figure |l| decrease with increasing m and have the limiting 
values 



lim g' n 



(m) 



(2.17) 



Zajc has concluded on the basis of this behaviour that || 
"the two-particle correlation function becomes a super- 
position of terms with successively broader distribution 
in Pi — P2, leading to an increasingly smaller value for 
the inferred radius." The reason is that Bose-Einstein 
symmetrization effects enhance the probability of finding 
bosons closer together in configuration space and hence 
result in broader q-distributions of the two-particle spec- 
tra. 

The above arguments explain the effect of multiparti- 
cle correlations qualitatively For a quantitative under- 
standing, the weights of higher order terms contributing 
to the one- and two-particle spectra are important. These 
weights are governed by the terms C m which (up to a cor- 
rection factor (1 — h)/ det B m of order unity) essentially 
decrease like (to — l)th powers of the inverse phase space 
volume. For event multiplicities in the hundreds, a quan- 
titative analysis can th en be done nume rically, using the 
analytical expressions ( 2.14 ) and ( |2.15|) . We defer such 
a study to a slightly more general model in section |^ 
where certain analytical simplifications allow for a more 
transparent discussion. A short comparison of the qual- 
itative and quantitative properties of the models stud- 
ied here and in section ^ is given in the text following 
Eq. @. 



III. GAUSSIAN WAVEPACKETS 



In the example of section II B, we have started from 

whose momentum 



single particle creation operators 



4 



support <?(k) does not carry a label i. As a consequence, 
the iV-particle states considered in section [I B are build 



up from N single particle wavefunctions with identical 
phase space localization. We now adopt a more general 
setting in which a set of N phase space points (pj, rj, ti) 
is associated with N Gaussian wavepackets centered 
at initial time U around the points (pj,rj) with spatial 

width a, HEHni 



fi(*,U) 



1 



(71-a 2 ) 3 / 4 



exp 



-^(X-r^+ip^X 



(3.1) 

is proportional to 
cxp [— cr 2 (k — pi) 2 /2 — in ■ (k — p.^] and can be com- 



The Fourier transform of /, 
-ir. 

pared to the momentum support of $ in ( [2. 11a ). The 
corresponding TV-particle symmetrized states reduce to 
those considered in section II B for po = 1 ja and pi — 
= rj, when the momentum support function becomes 
independent of the particle label i. 



The one boson state (3.1) is optimally localized around 
(pi,Ti) in the sense that it saturates the Heisenberg un- 
certainty relation Ax • Ap x = 1 , with Axi = a at initial 
time t = ti. Different choices of single particle wavefunc- 
tions can lead to different results for the calculated spec- 
tra. We have some control over the extent to which differ- 
ent choices matter by finally checking the u-dependence 
of our results. The model parameter can range between 
a G [0,oo]. In Refs. jnUTBI], it was argued that a realis- 
tic value for a is for pions of the order of the Compton 
wavelength. 

The free time evolution of /j is specified by the one- 
particle hamiltonian Hq which acts as a multiplication 
operator in momentum space, 



(27T) 



^3 J d 3 ke ik.X- iEk ( t - ti )J. (k) 



(3.2) 



The resulting building blocks for the iV-particle spectra 
are 



iEpti 



A(P) = 2 3 / 2 (7ra 2 ) 3 / 4 e-( ff2 / 2 )(P*- p ) e 
«,(P) = ^i(P) = 2 3 (trt 2 ) 3 / 2 e- 2 ^- p ) 



(3.3a) 
(3.3b) 



To streamline our notation, we have neglected in Di a 
factor exp[ir^ • p,]. The product of these fa ctors cancels 
in the calculation of the Wigner function ( |2.1b| ) and a 
fortiori in all the functions derived from it. The function 
Sj(P) = Tu{P) denotes the probability that a boson in 
the state fi is detected with momentum P. This mea- 
sured momentum P has a Gaussian distribution around 
the central momentum p^ of the wavepacket. 

The functions in ( 2.4c ) characterize the overlap be- 
tween the wavepackets fi and fj and play an important 



role in the ZP-algorithm. They take a particularly simple 
form if all particles are emitted in a flash, 



fij oc exp 



x exp 



(3.4) 



All terms contributing to V]^ or contain the same 
number of factors fij and hence, the normalization of /y 
does not matter in what follows. For notational conve- 
nience, we change it by fixing fa = 1. The functions fij 
measure the distance \z% — Zj\ between the phase space 
points i and j. This distance measure depends on the 
wavepacket width a but leaves the phase space volume 
independent of a, 



\fij\ = exp [-\\zi - Zj\ 2 ] , 
1 



(3.5) 



IV. MULTIPARTICLE CORRELATIONS FOR A 
GAUSSIAN MODEL 

We now determine quantitatively multiparticle corre- 
lation effects for a source of N identical bosons whose 
wavepackets of spatial width a are emitted instanta- 
neously according to a Gaussian iV-particle phase space 
distribution (|2.6|) with 



S(t) 



7T 3 i? 3 A 3 



exp 



r 

fi 2 



P 

A 2 



(4.1) 



Our main aim is to study for this model the extent to 
which multiparticle correlations modify the slope of the 
one-particle spectrum and the width of the two-particle 
correlator. To this end, we calculate first the building 
blocks of the ZP-algorithm. Having explicit expressions 
for <7q™' , gfp and C m in terms of simple polynomials will 
simplify our discussion considerably. The corresponding 
first order Pratt terms are 



Gi(Pi,P 2 ) 



1 , 
(1 



(4.2a) 



a 2 A 2 ) 



-3/2 



x exp 



x exp 



R 2 



T T ' q 



-K 



A 2 + l/a- 



(4.2b) 



All higher order Pratt terms can be calculated explic- 
itly as averages over relative and average pair distribu- 
tions. Details are given in Appendix |b| The momentum- 
independent terms read 



5 



(^ m) 4 m) )- 3/2 



[i 



-][i 



2a 1 



-3(m-l)/2 



(m) _ 



1 

E 

fc=0 



771- 



(m) 



E 

fc=0 



1 



l + 2<7 2 /i? 2 ' 



1 + 2/<t 2 A 2 



(4.3a) 



(4.3b) 



(4.3c) 



(4.3d) 



Here, denotes the greatest integer not larger than | 
('floor'), and the least integer not smaller than | 
('ceiling'). The notational shorthands a and 6 range be- 
tween and 1 depending on the phase space localization 
of the wavepacket centers, i.e., they map the whole pa- 
rameter space 0<i?<oo,0<A<ooof the model ( |4.lD 
onto a finite region. The momentum-dependent terms are 
given by 



G m (Pi, P2 



2tt 
A 2 



(m) 



b 9 K 



x exp 



x exp 



(m) 

9Q 

(rn) 

9k 



hi m) /h. 



V4' 



(m) 



l + ]T(a&) fe 



fe=i 



+ 1 

(m) 



rn 
2fc 



3/2 



2) q 2 



(4.4a) 
(4.4b) 
(4.4c) 
(4.4d) 



The comparison of the present model calculation with the 
Zajc model ( 2.12 ) is not straightforward. As mentioned 
in the sequel of Eq. (3.1), the wavepackets used in both 
models can be compared by setting the wavepacket cen- 
terspi = = r,-. However, the integral over p(x) in 
( [2.11b ) performs an average over the positions x while 
in the model (4.1) we do not average over the position 



X but over the centers of wavepackets pi, rj. Due to 
these different starting points, the Zajc model is not a 
simple limiting case of (4.1). Nevertheless, main features 
of the Zajc model (2.12) can be reproduced qualitatively 
in the present model. The leading contribution of the 
momentum-independent Pratt terms C m shows again a 
power law behaviour. Also, the m-dependence of the 
terms ffg™ and g^ 1 of the Zajc model is recovered in 
certain limiting cases, 



lim lim o>, 
lim lim q^ 1 



(m) _ 



(4.5a) 
(4.5b) 

In general, however, the m-dependence of the terms g^ 
and g# is much weaker in the present model. Especially, 



there is no limit in which both g^ = ^ and g^ 1 ' = m. 
These differences between both models may provide a 
first idea about the extent to which the choice of the 
model distribution affects our conclusions. 



(m) 



A. Weighting multiparticle contributions 

The normalization u(N) is not a direct physical ob- 
servable, but it determines the weights with which mul- 
tiparticle correlations contribute to the particle spectra. 
To see this, we consider the one-particle spectrum, 



N 



"mG m (P, P)/ C„ 



(4.6) 



The m-th order contributions G m /C m are normalized to 
one, and the weights v m add up to unity, 



JY 



i = E 



_ (N — 1)1 uj{N-m) 
Vm -(N-m)\ lu(N) °" 



(4.7a) 
(4.7b) 



The lowest order contribution G\ /C\ , by which the one- 
particle spectrum is typically approximated, contributes 
a fraction v\ only, the value (1 — v±) characterizes the 
importance of higher order contributions. For a quanti- 
tative analysis we now determine the dependence of the 
normalization u>(N) and the weights v m on the event mul- 
tiplicity N and the phase space density of the emission 
region. 

We consider the terms C m , the building blo cks o f 
u>(N). For the present model, these are given in ( 4.3a ). 



The factor (h^h^) in this equation ranges between 1 
and 2 2 (" l ~ 1 ), and can be written as 



J corr ^ 




2(m-l) 



(4.8a) 
(4.8b) 



The correction factor fcml appears only linearly in the 
expressions for C m , rather than as an m-th power, and 
it is of order O(l) (it is exactly /i™r = 1 for a choice of 
parameters R and A such that a = b). This allows for 
the approximation 



C m - 
e(i?,A) = ([l 



ab\ 



-3/2 



(4.9a) 
(4.9b) 



with which the normalization uj(N) takes the simple form 



G 



JV 



k=l 



(fc) JN-k) 



N 



N 



r(-i + i) 



ri( i + e ( fc - 1 )) • 



(4.10) 



fc=i 



(k) 

Here, the combinatorical factors denote the number 
of permutations of N elements which contain exactly k 
cycles. They are commonly refered to as Stirling numbers 
of the first kind. We have used their generating function 
in terms of T-functions |25| . 

We can now determine the weights of multiparticle con- 
tributions to the one- and two-particle spectra. For the 



one-particle spectrum, we find by inserting (4.10) into 
the ZP-algorithm 



Pvoi 
Vi 



N e. 



1 



l + e(JV-l) 

Pvoi 
(1 + Pvoi)™ ' 



Pvoi 



(4.11a) 
(4.11b) 

(4.11c) 



The approximation in the last line is valid for large mul- 
tiplicities, when m <C N . Similarly, the weights u m for 
the different contributions to the two-particle momen- 
tum spectrum can be calculated. Using the power law 
behaviour C m — e m_1 , wc find 

JV TO— 1 

7^(Pi,P 2 ) =J2 Um J2 ffi >m -i(Pi,P 2 ), (4.12a) 

rn—2 i—l 

Gi(Pl,Pl) G m _i(P 2 ,P 2 ) 



Hi, m -i(Pi, P2) 



Ci 



G l (Pi,P 2 )G m ^(P 2 ,P 1 ) 



a 



G„ 



(4.12b) 



where 



(JV-2)! w(N-m) 



(N-m)\ u(N) 

V m -1 



TO — 2 

Pvol 



(l + Pvol) m (1 + pvol) 



(4.13) 



Again, the approximation in the last line is valid for large 
multiplicities, when m <C N . To sum up: multiparticle 
correlations account for a fraction p vo i / ( 1 + Pvoi) of the 
one-particle spectrum. For the two-particle spectrum, 
they are somewhat more important: the pure pairwise 
correlations receive only a weight 1/(1 + p vo i) 2 - 

For sufficiently large event multiplicities, the weights 
V m and u m of multiparticle contributions are not sepa- 
rate functions of e and N, but functions of the product 
Pvoi only. The physics entering p vo \ can be most easily 
illustrated in the large phase space volume limit, when 



(i? 3 A 3 ) ' 



R 

for — , A • a > 1 
a 



(4.14) 



We hence call p vo \ a "phase space density of emission 
points" . This notion should not be taken too literally: 
the product of the volumes of three-dimensional spheres 
in position and momentum space is (|-7r) 2 i? 3 A 3 , and 
hence, p vo \ is for large sources approximately a factor 
10 larger than the particle number per unit phase space 
cell. Also, for realistic source sizes, the value of e devi- 



ates significantly from the approximation (4.14), and a 
calculation of p vo \ starting from (4.9b) is preferable. 

One can ask whether in the large N limit, the nor- 
malization ui(N) becomes a function of p vo \ only. To 
clarify this, we recall that a product Y\k=i 0- + a fc) w ^ tn 
cik > has a N — > oo-limit if and only if X^fcli °fe con- 
verges. For the normalization ( 4.10 ), we find Y2k=i ° fe 
= I cN (N — 1), i.e., for fixed phase space density, 



N 



lim uj(N)\ f 



fixed 



(4.15) 



There is no physical reason why lo(N) should remain fi- 
nite. It is not an observable. What matters in a quan- 
titative study of multiparticle correlation effects is the 
phase space density of emission points and not the par- 
ticle multiplicity, what matters are the weights v m and 
u m , and not the normalization lo(N). 

We finally estimate realistic values for the phase space 
density p vo i in heavy ion collisions. For a choice of model 
parameters R « 5 fm, a w 1 fm, A ps 150 MeV say, 
we find e ~ 10~ 2 . With multiplicities of like-sign pions 
in the hundreds, this leads to p vo i > 1. The present 
static, spherically symmetric model is however unrealistic 
in so far that it does not take the strong longitudinal 
expansion into account which significantly increases the 
volume out of which particles are emitted. From these 
heuristic considerations, we expect realistic phase space 
densities to lie in the range 



0.1 < p vol = eN < 1.0. 



(4.16) 



Depending on the precise value of p V oi in this range, the 
importance of multiparticle contributions to the one- and 
two-particle spectrum varies significantly. For p vo \ = 1.0, 
higher order contributions start dominating, while they 
account for sa 10 % of the signal if p vo i =0.1. 



B. 



The momentum dependence of multiparticle 
contributions 



The two m-dependent terms 
trol the dependence of G m (Pi,P 2 ) on the relative pair 
momentum q and the average pair momentum K. They 
determine the momentum dependence of all particle spec- 
tra. In Figure |^, these factors are shown as a function 
of b for a source with 7? = 5 fm and a = 1 fm. The 



and Qq^ in (4.4) con- 
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first message of this plot is that even for very high mul- 
tiparticle contributions (e.g. m = 100), the momentum 
dependence of all building blocks G m can be calculated 
exactly. Secondly, the factors gg and bg^"' show the 
interesting property that irrespective of the value cho- 
sen for A (and hence for 6), they rapidly converge to 
an m-indcpcndent quantity. In contrast to the limiting 
case ( 4.5a ), the m-dependence of <7q is much weaker for 



realistic model parameters a, 6, 



9Q 



(4.17) 



Analogously, for realistic model parameters a, b, 

bg^-^bgK, (4.18) 



while for the limit (4.5b) of the parameter space, a strong 
m-dependence remains. 



convergence to the limiting behaviour (4.17), (4.18) how- 
ever is observed for all values a > 0.1. For choices of the 
model parameters realistic for heavy ion collisions, the 
factors b <?#™ and <?q are hence well approximated by 
an m-independent constant for sufficiently large m. We 
can write the higher order Pratt terms as 



G m (Pi,P 2 ) = C m f g exp[-Aq 2 - BK 2 ] , (4.19a) 



/, = ( %bg K 



3/2 



A = 



r 2 R 2 



9Q 



R 2b 

B=^9K- 



(4.19b) 
(4.19c) 
(4.19d) 



In subsection IV D| , we exploit this simple m-dependence 



of the terms g^\ g^y 1 



0.2 0.4 0.6 0.8 1.0 




0.2 0.4 0.6 0.8 1.0 

FIG. 2. Numerical calculation of the factors g^ andbg^ 
which determine the q- and K- dependence of the m-th order 
Pratt terms. For a source with spatial radius R — 5 fm and 
a wavepacket width a = 1 fm, the plot shows and b fljf 
as a function of b — 1/(1 + 2/<j 2 A 2 ). Different lines denote 
different orders, m — 2 (dash- dotted), m = 3 (dotted), m — 5 
(dashed), m = 10 (thin solid) and m — 100 (solid). 



We have checked that the conclusions ( 4.17 ), ( 4.18| ) 
hold for a large range of the model parameter space, in- 
cluding the part realistic for heavy ion collisions. Going 
to smaller source sizes R (and hence to smaller values for 
a), the factor gQ is found to deviate significantly from 
unity. Also, values for bgn vary significantly. The rapid 



C. The one-particle spectrum 

For the discussion of the one-particle spectrum, we in- 
troduce the temperature T via 



2MT . 



(4.20) 



The model p(p,r) in (4.1) describes then a phase space 
distribution of emission points with Boltzmann tempera- 
ture T . Our aim is to determine how the slope and shape 
of the observed one-particle spectrum V N (P) changes 
with the slope T of this distribution p(p, r) and to what 

extent it is affected by multiparticle correlations. ■pjv(P) 
is a superposition of Gaussians of different widths 



^(P)oc^ Wm e-^/ T -(™), 

1 



N 



T eS {l) = T^ n ' = T 



Teff(m) = 



T 



2b g^ 



2Ma 2 ' 
m > 1 , 



(4.21a) 

(4.21b) 
(4.21c) 



where E-p — P 2 /2M, and T G ff (m) characterizes the slope 
of the m-th order contribution. According to (4.21), the 
one-particle spectrum cannot be characterized by a single 
slope parameter. For a qualitative understanding, we 



consider first the largest slope parameter T P g lr = T c ff(l) 
and the smallest slope parameter TJg ult = T e gf(m ^> 1). 
Here, the superscripts pair and mult stand for 'pairwise' 
and 'multiparticle' correlations. 

If all multiparticle contributions vanish, then the mo- 
mentum distribution of V N (P) coincides with that of 
Gi(P,P). The slope T o p ff air of G*i(P,P) is shown in Fig- 
ure 0a as a function of the model temperature T for the 
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pion mass M = 139 MeV and different wavepacket lo- 
calizations a. T^g" is always larger than the model tem- 
perature T. For a spatial wavepacket width a = 1 fm 



e.g., the term 



takes a value of 140 MeV. Even for 



2M<r 2 

small model temperatures T as input, the quantum con- 
tribution 2Ma 1 accoun ts f° r a slope parameter T c ^ lr com- 
parable to the Hagedorn temperature. For this appar- 
ently leading effect, the notion 'quantum temperature 
was coined H. 



> 

CD 
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'a 3= 
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0- 



(a) 




...- = 0.7 fm 

a = 1.0 fm 
a = 2.0 fm 



50 
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150 200 

T in MeV 



o = 0.5 fm 




100 150 200 
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FIG. 3. The one-particle slope parameters T^g T and T^ lt 
characterize the limiting cases of vanishing and dominant 
multiparticle correlation effects. They are shown as func- 
tions of the model temperature T for different values of the 
wavepacket width o . The diagonals for T c g = T (thin solid 
lines) are included to guide the eye. 

Figure ||b shows that multiparticle contributions can 
change this picture qualitatively, if they are dominant. 
The slope parameter T™ 1 ' still depends significantly on 
the choice of the wavepacket localization. But for model 
temperatures in the range 100 MeV < T < 200 MeV, 
there is always a value for a, 0.7 fm < a < 1 fm, such that 
the observed temperature T^ ult coincides with the model 
temperature T. For sufficiently large a, the multiparticle 
effect can even overcompensate the broadening due to 
the quantum mechanical localization, TjS < T. This 
illustrates that multiparticle symmetrization effects tend 
to populate the low momentum region of the one-particle 
spectrum, thereby increasing the slope of the spectrum. 
Both effects, this narrowing and the broadening due to 
the quantum mechanical localization, are governed by 
the same scale a and hence they cancel at least to some 
extent. 

The lowest order term G\ with slope parameter T^g 1T 



contributes a fraction 1 / (l+p vo \) to the one-particle spec- 
trum only. Hence, the one-particle spectrum V N is not 
monoexponential, but can be be characterized by local 
slope parameters which specify the tangent onto V N at 
Ep = P 2 /2M. Going to very large values of Ep, the 
local slope always coincides with that of T^ n , since the 



exponential of slowest decrease dominates in ( |4.21a[ ) . For 
sufficiently small Ep (below 2 GeV say), however, and 
for realistic phase space densities, neither the pair nor 
the multiparticle contributions can be neglected. The 
local slope lies between T e p ff air and T c ^ ult . For a more 
quantitative statement, we have plotted in Figure ^ the 
one-particle spectrum for different phase space densities 
function of Ep. Increasing p vo \ and hence the 
contribution of multiparticle correlations, the local slope 
of the spectrum becomes steeper. Also, the superposition 
of Gaussians of different width leads to a slight curvature 
of the spectrum. Fitting the spectra in the range 0-1 
GeV naively with a monoexponential thermal distribu- 
tion, a variation of p vo \ between and 0.5 results in a 
change of the slope parameter 



AT, 



eff 



20 MeV. 



(4.22) 



This has to be taken properly into account in a quanti- 
tative analysis of hadron spectra. 



w 10 



800 




FIG. 4. The one-particle spectrum of a source typical for 
heavy ion collisions (R = 5 fm, T = 100 MeV, a — 1.2 fm) 
becomes steeper with increasing phase space density p vo \ of 
emission points. The solid line characterizes a monoexponen- 
tial behaviour and is included for comparison. 



We finally compare the spectrum ( 4.21a ) to that of 
a Bosc-Einstein distribution as obtained from an arbi- 
trary initial phase space distribution p(p, r), keeping the 
bosons in a box till they have equilibrated, 



1 



e (PV2M- M )/T _ l 



oo 

E 

TCI— 1 



m -m(P 2 /2M- /J )/T 
V BE C 



(4.23a) 
(4.23b) 



For the spectrum V N in (4.21) being of Bose-Einstein 
form, the weights v m and the Pratt terms G m would 
have to show a particular m-dependence, 
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const ■ v„ 
G m (P,P)/C v 



^BE 

exp[ 



;E P /T] 



(4.24a) 
(4.24b) 



In the present mode l, (|4.24a ) is satisfied by setting vbe = 
Pvoi/(l + Pvoi), see ( Hlb| ). However, the m-dependence 
of the terms G m (P, P) is in general weaker than what is 
required to match ( 4.24b| ). For the Zajc model, the in- 
dependence of G m is compatible with the Bose-Einstein 
distribution in the infinite volume limit. However, the 
Pratt terms C' m differ significantly from a power law 



which indicates a deviation from (4.24a). The reason for 

these differences is, that our model calculations of V ' N 
take a particular distribution p(p, r) as initial condition 
and include Bose-Einstein symmetrization, but they do 
not contain an equilibration mechanism: the particles are 
emitted and propagate freely. Also, the TV-particle state 
^ n is not an equilibriu m sta te. Hence, the spectrum V N 
depends in contrast to ( 4.23 ) on p(p,r). In general, it is 
not a Bose-Einstein distribution. 



j e -W-BK' 



N 



u(N) 

(IV - 2)! 
(N-J)l 



j-i 



u(N -J)J2Ci Cj-i , (4.27) 



J=2 

and the correlator reads 



c(Pi ,p 2) = 5W( 1 .- 2A '- 2BK 



u)(N) 



— BP 2 —BP 2 



(4.28) 



The normalization 



i(N) 



obtained in this approximation 
remains unchanged if the full m-dependence of the Pratt 
terms G m is included, though the momentum dependence 



of (4.28) is then much more involved. 



The structure of the normalization 



is important 



for what follows: The term u)(N) contains exactly Nl/2 
terms, while oj(N) contains TV! terms. By integrating 
u(N)Pff(Pxi P2) over P 2 and using the power law be- 



N 

haviour C. 



we find 



D. The two- particle correlator 

In this section, we first discuss the normalization of the 
two-particle correlator C(K, q) and how it is calculated 
in the so-called pair approximation. Then we turn to the 
study of multiparticle effects. 

1. Normalization of the two-particle correlator 

There has been some debate recently about the appro- 
priate normalization of the two-particle correlation func- 
tion used in the HBT analysis of multiparticle produc- 
tion. For a compilation of the different normalizations 
used, and their problems, see In the present work, 
we use the normalization Q 



TV 



C(Pi,P 2 



<7 n d 6 a 7rn /d 3 Pid 3 P 2 



(TV(TV - 1)) (d 3 ^/d 3 Pi) {d 3 a^/d 3 P 2 ) ' 

(4.25) 



where a„ is the total pion cross section. This origi- 
nates from normalizing both the single- and the double- 
differential cross sections to unity. Setting these cross 

sections proportional to the one- and two-particle spectra 
1 2 

V N (Pi), , P A r(Pi,P2), the two-particle correlator reads 



c(Pi,p 2 ) = _ffi P -i P2) 



(4.26) 



The simplest way to make further progress is to assume 
that the approximation (4.19) for G m is valid for all m > 
1. Then, the two-particle spectrum can be written as 



co(N)-Q(N) = Y, 



(N-2)\ 



j-i 



,7=2 



u>(N) 



(TV- J)! 
= 1 - e 



lo(N- 



J) Cj , (4.29a) 
(4.29b) 



i=i 



The normalization of the correlator is smaller than unity, 
the offset depends on the phase space volume occupied 
by the source. We shortly comment on the importance 
of this result: 

The two-particle correlator can be viewed as the factor 
relating two-particle differential cross sections obe of the 
real world (where Bose-Einstein symmetrization exists) 
to an idealized world in which Bose-Einstein correlations 
are absent, 



<tbe(K, q) = C(K, q)cr N o(K,q). 



(4.30) 



If C(K,q) > 1 everywhere, this implies that Bose- 
Einstein symmetrization effects increase the total cross 
section. For heavy ion collisions, there is no direct test of 
whether such an enhancement exists. In e + e~-collisions 
however, we know that Bose-Einstein correlations do not 
affect total cross sections appreciably, since e.g. pertur- 
bative QCD predicts the production characteristics of Zq 
to a per mille level without invoking them. In our calcu- 
lation, we start from an TV-particle state and we require 
that after final state Bose-Einstein symmetrization, N 
particles are detected. Hence, we explicitly assume that 
final state Bose-Einstein correlations do not affect the 
total cross sections. According to the above calculation, 
this automatically implies an offset of the normalization 



of C(K,q) below unity. The offset e in (4.29t) measures 
the system size from which particles are emitted. This is 
intrinsically consistent: for smaller system sizes (e large), 
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the correlator shows an enhancement in a broader q mo- 
mentum region, and the normalization (1 — e) is smaller. 



This ensures that Og^; 



'NO- 



We finally note that irrespective of the offset 
li(N) /w(N) of the normalization, the correlator in ( 4.28Q 
changes by a factor 2 between the limits q = and 
q — * oo. This is a consequence of the approximation 
( 4.19| ) which we have justified for the present model in 
subsection IV B. In general, the K- and q- dependence 
of the Pratt terms G m leads to a more complicated de- 



pendence of the correlator C(K,q) according to (4.12) 



Depending on the model, this may affect the intercept 
of C(K,q) at q = 0. Indeed, a decrease of the inter- 
cept with increasing event multiplicities was reported in 
recent model studies |§,|lo|. In contrast to the present 
study, these models however do not work with fixed event 
multiplicities and show a significantly different physics, 
including pion lasing effects [Q. We only conclude from 
the present study that multiparticle symmetrization ef- 
fect do not lead automatically to a strong decrease of the 
intercept parameter. 



2. The normalization in the pair approximation 

We now explain why the offset (1 — e) of the normal- 
ization of C(K,q) is not obtained in conventional cal- 
culations where multiparticle symmetrization effects are 
neglected. For the wavepackets introduced in section [II, 
the wavepacket overlap functions fij 



in (3.4) shows a 
Gaussian decrease with the phase space distance \z% — Zj\. 
In general, the 



pair approximation: 



(4.31) 



is an uncontrolled approximation to this Gaussian be- 
haviour. It becomes exact in certain limiting cases 



lim fij = lim /, 



■ i — 



(4.32) 



In the pair approximation, all higher order Pratt terms 
vanish 



C m = G m (K,q) =0, m>2. 
for fij =$ij, 



(4.33) 



and one regains the results of the conventional calcula- 
tions where multiparticle effects are neglected: the two- 
particle correlator is normalized to unity, 



to(N) = u(N) = 1 for f i: 



and ( |2.6| ) does not change the phase 
tribution since u(N) = ("J n) ■ 
Gi(Pi,P 2 )Gi(Pi,P 2 ) in f $I% th en insures 
C(Pi,P 2 ) > 1. According to (1§) 



(4.34) 

space dis- 
Thc term 
that 

this contradicts 



the statement 



eri?<~). The origin of this difficulty 



can be traced back to the integral 



d 3 Pxd 3 P 2 Gi(Pi,P 2 ) 



G 2 



(4.35) 



which should vanish according to ( 4.33 ). The pair ap- 
proximati on do es not treat the integral G 2 and the inte- 
grand of ( 4.35 ) on an equal footing. It sets C 2 = but 
uses Gi(Pi,P 2 ) 2 for the calculation of the Bose-Einstein 
enhancement. It is this inconsistency which leads to a 
correlator G(K, q) > 1 everywhere. 



3. The HBT-radius parameters 

Once the momentum-dependent higher order Pratt 
terms G m and their weights v m , u m are known, the deter- 



mination of the two-particle correlator (4.26) is a matter 
of straightforward numerical calculation. From G(K, q), 
the HBT radius parameters are obtained by fitting the 
Gaussian ansatz 



G(K, q) oc 1 + Aexp [-i? H BTq 2 ] 



(4.36) 



We restrict our discussion to this onc-dimcnsional 
parametrization of G(K,q) since our model (4.1) is 
spherically symmetric. The general, analytical form of 
the t wo-pa rticle correlator ( 1.26] ) is obtained from (4.6) 
and ( 4.12 ) and it is quite involved. For a transparent 
discussion, we hence turn first to limiting cases. In the 
pair approximation, whe n all multiparticle effects are ne- 
glected, we obtain from ( 4.2b ) the well-known result |ll|] 



• 2 R 2 

□pair z _ ll 
"HBT — 2 



AV 2 



2 1 



(4.37) 



The other, equally unrealistic limiting case is that multi- 
particle contributions dominate completely. The result- 
ing two-particle correlator is given in (4.28) and the cor- 
responding HBT-radius parameter reads 



omult z 
n HBT 



^ a 2 g Q o- 2 A 2 + 2g Q - g K 
a 2 A 2 + 2 



-9Q 



(4.38) 



The difference between these two expressions is signifi- 
cant. For a discussion of the main qualitative and quan- 
titative effects, we now focus on the parameter region 
R 2 ^S> o 2 relevant for relativistic heavy ion collisions. 



I n th is regime, the ^--dependent parts of ( 4.37 ) and 
(4.38) can be neglected, and according to the study of 
section IV B, gn w 1. We then find the simple relation 



R 



mult 
HBT 



f ppair 
^2 HBT ) 



for 



(4.39) 



This is the factor 2, obtained first by Zajc |j) in his com- 
parison of the first and second order Pratt terms, see 
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given an estimate that in his model, 
multiparticle effects change the HBT radius parameter 
by as much as a factor 0.67 for a particle density of 1 per 
unit phase space cell ||. This would affect estimates of 
the source volume or the energy density by a factor 3. In 
Zajc's model, however, higher order multiparticle effects 
(to > 2) change the q-dependence of G(K, q) much more 
dramatically than in the present study, where has a 
very weak TO-dependence, negligible for R 2 3> cr 2 . 

In our model, the correlator is a weighted superposi- 
tion of Gaussians in q whose widths differ by as much 
as a factor \/2. In Figure 0, we have plotted the result- 
ing two-particle correlator for relatively low phase space 
densities where the lowest order Pratt term Gi is still the 
leading contribution. The result of the pair approxima- 
tion describes the main behaviour, but deviations due to 
multiparticle contributions affect the HBT-radius param- 
eter on a 10 %-level for moderate phase space densities. 
This 10% effect translates into underestimating the vol- 
ume by 30 % and overestimating the corresponding en- 
ergy density by a similar amount. In models which show 
a stronger TO-dependence of the momentum-dependent 
part of G m , the effect may be significantly stronger. In 
general, the degree to which multiparticle effects mod- 
ify HBT radius parameters depends significantly on the 
phase space density of emission points. 
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FIG. 5. Multiparticle symmetrization effects broaden the 
two-particle correlator. For a source size of R — 5 fm, the 
plot shows the resulting two-particle correlator for different 
particle phase space densities p voi . 



V. MULTIPARTICLE CORRELATIONS FOR 
ARBITRARY MODELS 

In this section, we shortly discuss how the calculation 
of multiparticle correlation effects can be extended to 
more realistic source distributions where the analytical 
techniques used in section IV are not applicable. Pratt 
has shown already how to calculate the terms C m and 
G m (Pi,P 2 ) for continuous source functions Q: up to 



be used, and an improved Metropolis method allows to 
push numerical calculations up to m « 20. This seems 
to be sufficient for all practical purposes, since according 
to our model study, the weights v m , u m for to > 20 are 
negligible for realistic phase space densities p vo i. 

Here, we propose yet a different technique which is 
applicable to events characterized by a set of N discrete 
phase space points (p^, r^ti). We interprete these points 
as the centers of Gaussian single-particle wavepackcts, 



i.e. 



(Pi,ri,ii) — > fi . 



(5.1) 



An arbitrary distribution of N discrete phase space 
points is the most general 'source model' in the present 
framework. We first discuss algorithms which allow for 
the calculation of the two-particle correlator from a set 
(pi,rj,ij). Then we comment on applications of these 
algorithms to event generators. 



A. An algorithm for Bose-Einstein weights 

We start with the simplest case, the pair approxima- 
tion, which reduces the sum over N\ terms in the two- 
particle spectrum to a sum over all ±N(N - 1) particle 
pairs (i, j). Each pair is weighted with the pair probabil- 
ity Vij, calculated from the corresponding two-particle 
symmetrized state, O] 



^(P lf P 2 ) = i 



D i (P 1 )Z> i (P 3 )+A(Pa)£ J -(Pi) 
2 



N(N - 1) 



n 5><,(Pi,Pa) 



(5.2) 



order m < 5, straightforward Monte Carlo methods can 



To compare this expression to the ZP-algorithm, we in- 
troduce the quantities 

G 1 (P 1 ,P 2 ) = i£ti £>*(Px)A(P 2 ) , (5.3a) 

T c (P lf P 2 ) = -^Eli A*(Pi)A(Pi) 

x D*(P 2 ) Di(P 2 ) . (5.3b) 

Here, Gi(Pi,P 2 ) is the first order Pratt term for the 
discrete phase space distribution (pi,rj,tj) and T c is a 
finite multiplicity correction, correcting for the double 
counting of identical pairs (i, i) in 

^P ij (P 1 ,P 2 ) = G 1 (P 1 ,P 2 )G 1 (P 1) P 2 ) 

+G 1 (P 1 ,P 2 )G 1 (P 1 ,P 2 ) 
-2T C (P 1 ,P 2 ). (5.4) 

In the large multiplicity limit, when T c can be neglected, 
this expression is equivalent to the pair approximation 



of the ZP-algorithm, see e.g. (4.12a). The one-particle 
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spectrum v and the 2-particle correlation obtained in this 
way have been discussed already extensively in jflj . They 
can be rewritten in terms of the single-particle probabil- 
ities Si of (3.3b), 



C(K,q) = l 



e — V/2 



EL *i(K)e*«-» 



i/(Pi)i/(Pa)-T c 



AT 



Kp) = ^>(p), 

1=1 

iV 

r c (Pi,p 2 ) = ^« i (Pi)« i (p 2 ) 



(5.5a) 
(5.5b) 



The number of numerical operations needed to calcu- 
late this correlator grows linearly with the multiplicity 
N rather than quadratic as in ( |5.2| ) and this makes it 
particularly suitable for a numerical algorithm. At given 
observed momenta K, q, one calculates for each event the 
one particle probabilities Si(K) according to ( |3.3b| ) and 
performes the sums in ( [5.5| ) . The s, (K) are continuous in 
the measured momentum, and hence the obtained corre- 
lator is a continuous function of K and q, i.e., no binning 
is required. At high multiplici ties , when the factors T c 
are negligible, the correlator (|]^) can be rewritten as 
a Fourier transform over an emission function S(x,1C), 
thereby regaining the well-known starting point of most 
model studies, (see Refs. p| , ^5| for further details) 



C(K,q) = l 



Jd 4 xS(x,K)e i: 



(5.6a) 



f&xSfaPi) Jd*yS(y,P 2 ) 

OO 

S(x,K)=J2 S i( x > K ), ( 5 - 6b ) 

i=l 

Si(x, K) oc S(t - U) e~M x - r ') 2 - a2 ( K - P^) 2 . 

The above discussion shows, that in the pair approxima- 
tion, a numerical algorithm for the calculation of one- 
and two-particle spectra can be based on a discrete ver- 
sion of the first order Pratt terms G\. Here, we propose 
to extend this approach to take multiparticle correlations 
into account by calculating 

_ (N-m)\ ^ 



N\ 



G m (Pi, P 2 ) — 



(AT -to)! 



E ££(Pi)/m 



i% ...i m 

■/wii m A m (P 2 ). 



(5.7b) 



The sums in these expressions run over all sets of to out of 
N integers and over all ml permutations of each set. This 
implies that for event multiplicities in the hundreds, only 
terms up to order m«5 can be calculated in a reasonable 
amount of CPU time. A tentative strategy for calculating 



multiparticle correlation effects on the basis of ( |5.7| ) then 
proceeds as follows: 

1. Fit the calculated terms C m to a power law C m — 



and use the parameter e thus determined for 
a calculation of the weights v m , u m in (fLq) and 



(4.12a) 



2. Fit Gaussians in K and q to the momentum- 
dependent Pratt terms G m (Pi,P2). This allows 
to extract the terms and which govern 
the momentum dependence of multiparticle effects. 

3. Calculate the one- and two-particle spectra by in- 
cluding up to the numerically determined order, 
m = 5 say, all contributions exactly, approximating 
the momentum dependence of higher order terms 



by setting g^ 



(5) (m) 

9k\ 9q 



,( 5 ) 



for to > 5. 



This scheme draws on the experience gained from the 
study of the Gaussian toy model in section We expect 
that it allows to estimate multiparticle contributions for 
model distributions (pj,rj,tj). 



B. Bose-Einstein weights for event generators 

Many event generators for the simulation of heavy ion 
collisions have been developed in recent years |lrj-[20[|, 
and irrespective of the large variety of physical inputs 
present in their codes, the typical output of their event 
simulation contains for each event a set S of N phase 
space points (pj,rj,tj) which one associates to the final 
state particles produced. However, a choice of interpreta- 
tion is involved in comparing an event generator output 
z i = (Pii r ii^i) to experimental data. Usually, the mea- 
sured one-particle spectra are compared directly to the 
binned momenta p^, i.e. one implicitly interprets the 
simulated phase space points Zi as defining momentum 
eigenstates, which a fortiori carry no space-time infor- 
mation. For the calculation of two-particle correlations, 
this interpretation is not suitable, it leads to a sharp 6- 
like correlator jllj]. Also, a classical interpretation of 
Zi which takes both r.^ and p^ as "sharp" information, 
is problematic: ignoring the correct quantum mechani- 
cal localization leads to quantitatively and qualitatively 
unreliable results (^3|. On the other hand, including a 
quantum mechanical localization width tr, one changes 
both the one- and two-particle spectrum by a prescrip- 
tion which has no dynamical foundation. 

The origin of these difficulties in finding a consistent in- 
terpretation of the event generator output Zj = (pj, rj, ti) 
is well-known, see e.g. f2l"| , |22|| : quantum mechanical pro- 
cesses require a description in terms of amplitudes, while 
event simulations are formulated in a probabilistic set- 
ting. Bose-Einstein correlations occur by symmetrizing 
the production amplitudes of identical particles and are 
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hence not encoded for in event generators. The mere fact 
that the simulated output is a discrete phase space dis- 
tribution of emission points without identical multipar- 
ticle correlations (rather than a set of detected momenta 
which includes all multiparticle correlations and hence 
all space-time information available from the measure- 
ment), is the very consequence of using a probabilistic 
language. In practice, this forces one to take recourse to 
an algorithm which associates a posteriori Bose-Einstein 
weights with the simulated phase space distribution E, 
rather than obtaining these correlations from the dy- 
namical propagation of properly (anti)-symmetrized N- 
particle states. This a posteriori modification of an in- 
complete quantum-dynamical evolution creates the in- 
terpretational problems of the output z% = (pi,rj,tj). 
Phcnomenologically motivated numerical simulations of 
heavy ion collisions have remarkable success despite these 
fundamental problems. It is of some interest to confront 
them with experimental two-particle correlations, since 
this provides a test of their spatio-temporal (rather than 
only their momentum-dependent) prop erties . We expect 
that the algorithms discussed in section V A are useful in 
such comparisons. 



VI. CONCLUSION 

The Zajc-Pratt algorithm provides a simple technique 
for the calculation of multiparticle symmetrization ef- 
fects in one- and two-particle spectra. Based on this 
algorithm, we have studied to what extent multiparti- 
cle correlations steepen the slope of the one-particle and 
broaden the width of the two-particle spectrum. The 
scale of both effects depends sensitively on the particle 
phase space density p vo \ in the emission region. Multipar- 
ticle correlations contribute a fraction p vo i/(l + p vo \) to 
the one-particle spectrum and a fraction 1 — 1/(1 + p V o\) 2 
to the two-particle spectrum. Also, the m-dependence of 
higher order Pratt terms, i.e., the extent to which the mo- 
mentum dependence of multiparticle contributions G m 
changes with increasing order m, plays an important role. 
In a Gaussian source model with instantaneous emission 
and moderate particle phase space densities, the slope 
parameter T e g of the one-particle spectrum changes by 
up to 20 MeV, and the change in the HBT radius param- 
eters is of the order of 10 percent, if multiparticle effects 
are neglected. 

For the reconstruction program advocated in [jlj, our 
results indicate that neglecting multiparticle effects, one 
underestimates the source temperature significantly and 
one overestimates the energy density by up to 30 %. 
Here, the caveat is however that our calculation does not 
include multiparticle final state Coulomb interactions. 
Bose-Einstein and Coulomb effects typically arise on the 
same scale and compensate each other at least to some 
extent. This may significantly reduce the effect of multi- 



particle correlations in the measured ir + and tt~ spectra 
on which most of the phenomenological analysis is based 
currently. It will hence be very interesting to compare the 
slopes of the one-particle spectra of charged and neutral 
pions p?| |. On the basis of the above heuristic ideas, one 
may expect the slope of the 7r°-spectrum to be some- 
what steeper, since multiparticle Coulomb interactions 
cannot compensate for the multiparticle Bose-Einstein 
symmetrization effects. 

Our analysis of multiparticle effects has also shed new 
light on a recent discussion about the normalization of 
t he t wo-particle correlator. Defining the correlator via 
( 1.25| ), it is normalized for large relative momenta q to 
(1 — e), and not to unity. The i mportan ce of this result 
was already discussed in section [VDl: the normaliza- 



tion ( 4.29b ) corrects the result of the commonly used pair 
approximation and leads to Bose-Einstein unit weights 
which conserve event probabilities. 

We close by embedding the present work into a wider 
perspective: Many aspects of the dynamical evolution of 
relativistic heavy ion collisions are mesoscopic, and this 
makes it very difficult to decide whether certain physical 
observables have a known conventional interpretation or 
are indicative of the new in medium properties which the 
current experimental programs at CERN and RHIC aim 
for. The physics of final state hadrons however, though 
being mesoscopic and hence difficult to treat, is known 
in principle. It allows hence for a detailed quantitative 
description on the basis of present day knowledge. In 
a second step, such a description will pave the way for 
quantitative statements about the geometry and dynam- 
ics of the final stage of the collision process, thereby set- 
ting up an experimentally motivated starting point for 
a dynamical back extrapolation We hope that the 
present work, by quantifying one mesoscopic multiparti- 
cle effect, will prove useful in refining this reconstruction 
program and in separating known physics from the in 
medium properties we are looking for. 
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APPENDIX A: DERIVATION OF THE 
ZP- ALGORITHM 

So far, there is no self-contained derivation of the Zajc- 
Pratt algorithm in the literature. For the convenience of 
the reader, we give here the main arguments. 

Notation: We denote the term fi x i 2 by an arrow from 
i\ to %2- m of these terms form a closed TO-cycle C rn . The 
only other structure appearing in all spectra are open re- 
cycles Gmi see Figure |(| Here, D* is represented by 
a cross from which the arrow starts, by a cross at 
which the arrow ends. The Pratt terms C m , G m are 
obtained by averaging the emission points in Figure || 
over a model-distribution p(p, r, t). Dots without indices 
attached indicate that this average was carried out, see 
Figure fj]. 






G^,P 2 ) 

FIG. 6. Diagrammatic representation of the types of factors 
contributing to multiparticle spectra. 

— 

Normalization: The normalization w(N) — V N sums 
the products Cj 1 C l 2 2 ...C*» over all partitions (n, Z n )iV- 
Each partition of (n,Z„)jv has N\/ J7„ nl ln (l n \) possible 
realizations. These have to be multiplied by the (n — 
1)! different ways to combine each set of n points to a 
closed cycle. This leads to the prefactor TV!/ JT n ln (l n \) 
in (|~ 



(0(1)= • 

(0(2) = • • + 

ro(3) = • • • + 3 • m-»m 



FIG. 7. Diagrammatic representation of the averaged total 
multiplicity w(N). The 2-cycle stands for C2, the 3-cycle for 
Ci, etc. 

There exists a simple diagrammatic prescription to 
construct from the terms of w{N) the terms of w(N + 1). 
One adds to each diagram of w(N) an (N + l)th dot by 

• 1.) placing the dot once disconnected from all other 
cycles. 



• 2.) changing closed TO-cycles C m into to closed (m+ 
l)-cycles C m+ \. 

This prescription can be checked e.g. in Figure It 
automatically insures, that w(N) is represented by N\ 
diagrams, the maximal value for w(N) is N\. 

One-particle spectrum: Each term of u>(N) V N (P) has 
the structure w(N — m) G m : it contains exactly one open 
rn-cycle G m (P, P) , the N — m other dots being contained 
in closed cycles. To construct from it the terms of u(N + 
1) V N+1 (P), one has to use the replacement 



w(N - to) G n 



w(N + 1 - to) G m 
+w(N — m)m G m +i 



(Al) 



Diagrammatically, this is realized by employing the two 
diagrammatic rules given above and supplementing them 
with 

• 3.) change open TO-cycles G m into to open (m+ 1)- 
cycles G m +i- 

The prescription can be checked e.g. in Figure [§. Using 
the recursion relation ( |Al[ ), one can prove by complete in- 
duction that the averaged one-particle spectrum V N (P) 
takes the form (2.9). 



(o(D 

(0(2) 7T=* • + *Ht 

(0(3)^ 1 =)^« • f 

+ 2 * 



2* 



FIG. 8. Diagrammatic representation of the 1-particle spec- 
trum for multiplicities N = 1, 2, 3. 

2 

Two-particle spectrum: In w(iV) 7 3 Ar (Pi, P2), each 
term contains two open cycles of lengths i and j, the 
N — i— j remaining dots being contained in closed cycles. 
On the endpoints of the two open cycles Gi and Gj , the 
momenta Pi and P2 can be attached in different combi- 
nations: for i — j, there are two possiblities, for i =/= j, 
there are four. This combinatorics is properly encoded in 
the definition of the product Gi j = CiCjHij /(l + Sjj), 
see Eq. (4.12b). An N — > (N + 1) recursion relation for 
the two-particle spectrum is obtained by using the rule 
to change open TO-cycles into to open (to + l)-cycles, and 
properly treating the combinatorics for i = j, 



w(N -i-j) G itj 

— ► w(N -i-j) i(l + 5 i+ i t j - 5 it j)G i+ ij + j G lij+ i 
+w(N + l-i- j) G itj . (A2) 



The expressio n (|2.10| ) for the two-particle spectrum then 
follows from (|A2|) by complete induction. 
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APPENDIX B: CALCULATION OF 
MULTIPARTICLE EFFECTS 

Here, we give details of the calculation of higher order 
Pratt terms C m and G„, for the Gaussian source model 



(4.1) in section IV. Due to the structure of the functions 
fij , it is advantageous to change to relative and average 
momentum variables. For an average over to phase space 
points, we introduce the integration variables 



a,i=pi- p i+ i , h = ri- r i+ i , for i € [1, m - 1] , (Bla) 
a m = Vm + Pi , b m = r m + ri . (Bib) 

In terms of these variables, the building blocks for the 
m-th order Pratt terms read 



/l2/23"-/(m-l)rr, 
A r , 

frn ] 



e"-°i 



-Kb 2 



n 

3=1 

1 Y^ n— 1 I 1 v~^ m ~ 1 
~ 2 2^9 = 1 a 'j + 2 



1 b i A i , (B2a) 



j=n+l a 'J 



a 2 ( 

T (2^3=1 ' 



4tr 



(B2b) 



i ^-^.m— 1 , 



D 1 (P 1 )£>^(P a )=e 



x e 



xe * ^3 =1 *"■> , (B2c) 

-cr 2 (a m /2-K) 2 -iK^^" i 1 bj 
2 

1 ">> +iq ^. (B2d) 



a- 2 ( i v^m— 1 

--(q + E^i ■ 



For the Gaussian emission probability (4.1), the probabil- 
ity p2 of having two particles with wavepackets centered 
around phase space points (r^p^), (r,-,Pj), is given by 
the product p(r i; pj)- p(rj,Pj) and factorizes into a prob- 
ability distribution for the relative and average particle 
pair coordinates, 

p 2 (pi,pj,ri,r j ) = p(pi,T i )p(p j ,T j ) 

= Prel(Pi ~ Pj,r; - r,-,) 

Xpavc(P- i + Pj jTj + Tj) , (B3a) 

p(a, 6) = p rc i(a, 6) = p avc (a, b) 

1 _ a 2 _ t, 2 

e 2A^ is? . (B3b) 



(27r) 3 i? 3 A 3 



We then define the higher order Pratt terms as averages 
over pair distribution probabilities 

G TO (P!,P 2 ) = J ^jjtPajdPbjfta^bj) 

x£>l(Pl) /l2 /23 ... /(m-l)m A*„(P2) ■ (B4) 

This is a Gaussian integral which can be calculated ana- 
lytically. Its exponent is diagonal in all integration vari- 
ables bi and in a m , and doing the corresponding integrals 
leads to 



G m (Pi,P 2 )=([l + ^][l + £r] 



-3(m-l)/2 



-3/2 



q 2 e -i*(l-i)K 2 



x/(K,q), 

-l_(l + a6( m -l)) 



(B5a) 
(B5b) 



Here, the notational shorthands a, b of ( [4.3d ) are used 
again, 



7(K,q)=7r- 



„ / m — 1 
-3(m-l)/a y / TJ d 3 fl . 

x exp[ f),M,,« , + {wi + Vi)ri\ , 
iy s ; = q aVb , 

Vi = K 2<7V / 6a(2z — to) , 
My = % + b + a& (m - 1 - 2\i - j\ - 5 l0 ) 
a 2 b 2 



(B6a) 
(B6b) 
(B6c) 



(2i - m)(2j - to) . (B6d) 



a6(to - 1) + 1 
The integral I(K,q) can be calculated explicitly, 



dctM = 



/i 2 m) hj m) 
1 + ab(jn — 1) 



fuiM^/ujj =0, 



1 



( 7 - l)4 m )-6 7 ^ 

(m) ,(m) 



(B7a) 
(B7b) 
K 2 , (B7c) 



a 2 1 /i v " t; - W 
-WiMjy-Wj = - — , , 3 q 2 . 



(B7d) 



The factors , and /4"^ are defined in (4.3) and 
fi- 
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servables of model (4.1) can be shown [2] to depend only 
on the two effective source parameters i? 2 ff = R 2 + cr 2 /2 
and A 2 ff = A 2 + 1/2 ct 2 . The Pratt terms (4.3) and (4.4) 
on which our discussion is based, however, depend on 
all three parameters R, A and a. This shows that the 
calculation of higher order Pratt terms via pair distribu- 
tions (B4) destroys the one-particle factorizability. We 
withdraw the statement that the Pratt terms (4.3), (4.4) 
are solutions of the model distribution (4.1). They are 
solutions of (B4). 
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ERRATUM 

The results presented in section IV are not the results 
for the factorizable one-particle distribution (4.1), in con- 
trast to our statement. All results of section IV were cal- 
culated from the probability distribution (B3b) for rela- 
tive and average pair coordinates according to (B4). The 
use of this pair distribution (B3b) is equivalent to the 
one-particle distribution (4.1) only for the lowest order 
Pratt terms (4.2) on which the conventional two-particle 
symmetrized calculation of the Hanbury-Brown/Twiss 
effect is based. For the higher order Pratt terms, how- 
ever, our expressions (4.3), (4.4) differ from the results 
for the model (4.1), given in Ref. [1]. Especially, all ob- 
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